# File Information --------------------------------------------------------

  # Title: ReplicationAPublicDemandTheoryOfSanctions.R
  # Date: Spring 2025
  # R Version: R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
  # Purpose: Replicate the analysis for Caton and Webb (2025) - A Public
  # Demand Theory of Economic Sanctions, Manuscript and Appendix.

# Clear R -----------------------------------------------------------------
  rm(list = ls())

# Packages ----------------------------------------------------------------

  # Install Packages
  # install.packages(c("tidyverse", "xtable", "texreg","readstata13","margins","interplot"))

  # Load Packages
    library(tidyverse)
    library(xtable)
    library(texreg)
    library(readstata13)
    library(margins)
    library(interplot)

# Load Data ---------------------------------------------------------------
  
  # Working Directory
    setwd("Data/")

  # Load Data
    spbbdat <- as_tibble(read.dta13("DPDSspbb.dta"))
    glimpse(spbbdat)
    spbbdat
    
# Prepare Data ------------------------------------------------------------

  # Separate into Countries
    bb <- spbbdat |> filter(high == 0)
    sp <- spbbdat |> filter(high == 1)
    
# Rescale Political Variables ---------------------------------------------

  # Bolivia Brazil
    bb$handlingn <- bb$handlingn - 5
    bb$approvaln <- bb$approvaln - 5
    bb$competencen <- bb$competencen - 5
    
    sp$handlingn <- sp$handlingn - 5
    sp$approvaln <- sp$approvaln - 5
    sp$competencen <- sp$competencen - 5
    
    table(bb$handlingn) 
    table(bb$approvaln) 
    table(bb$competencen)
    
    table(sp$handlingn) 
    table(sp$approvaln) 
    table(sp$competencen)
    
# Figure 2: Competing Hypotheses about Supply and Demand for Sanct --------

  # Set Working Directory for Figures
    setwd("../Plots/")
    
  # Set Up Matrix
    Effect <- c("Demand Shift","Demand Elasticity")
    d <- c(20:1,20:1,sort(seq(5,14.5,.5),decreasing = TRUE),25:5)
    sc <- c(1:20,1:20,1:20,1:20)
    
    seq(5,14.5,.5)
    
  # Values for Figures
    hyp_dat <- as_tibble(cbind(d, sc, rbind(cbind(rep("Demand Elasticity",20)),cbind(rep("Demand Shift",20)),
                                            cbind(rep("Demand Elasticity",20)),cbind(rep("Demand Shift",20))),
                               rbind(cbind(rep("Before",40)),cbind(rep("After",40)))))
    hyp_dat <- hyp_dat |> rename("Effect" = V3)
    
    hyp_dat$Demand <- NA
    
    hyp_dat$Demand[hyp_dat$d == 1 ] <- 1 
    hyp_dat$Demand[hyp_dat$d == 2 ] <- 2 
    hyp_dat$Demand[hyp_dat$d == 3 ] <- 3 
    hyp_dat$Demand[hyp_dat$d == 4 ] <- 4 
    hyp_dat$Demand[hyp_dat$d == 5 ] <- 5 
    hyp_dat$Demand[hyp_dat$d == 6 ] <- 6 
    hyp_dat$Demand[hyp_dat$d == 7 ] <- 7 
    hyp_dat$Demand[hyp_dat$d == 8 ] <- 8 
    hyp_dat$Demand[hyp_dat$d == 9 ] <- 9 
    hyp_dat$Demand[hyp_dat$d == 10] <- 10
    hyp_dat$Demand[hyp_dat$d == 11] <- 11
    hyp_dat$Demand[hyp_dat$d == 12] <- 12
    hyp_dat$Demand[hyp_dat$d == 13] <- 13
    hyp_dat$Demand[hyp_dat$d == 14] <- 14
    hyp_dat$Demand[hyp_dat$d == 15] <- 15
    hyp_dat$Demand[hyp_dat$d == 16] <- 16
    hyp_dat$Demand[hyp_dat$d == 17] <- 17
    hyp_dat$Demand[hyp_dat$d == 18] <- 18
    hyp_dat$Demand[hyp_dat$d == 19] <- 19
    hyp_dat$Demand[hyp_dat$d == 20] <- 20
    hyp_dat$Demand[hyp_dat$d == 21] <- 21
    hyp_dat$Demand[hyp_dat$d == 22] <- 22
    hyp_dat$Demand[hyp_dat$d == 23] <- 23
    hyp_dat$Demand[hyp_dat$d == 24] <- 24
    hyp_dat$Demand[hyp_dat$d == 25] <- 25
    
    hyp_dat$Demand[hyp_dat$d == 26] <- 26
    hyp_dat$Demand[hyp_dat$d == 27] <- 27
    hyp_dat$Demand[hyp_dat$d == 28] <- 28
    hyp_dat$Demand[hyp_dat$d == 29] <- 29
    hyp_dat$Demand[hyp_dat$d == 30] <- 30
    hyp_dat$Demand[hyp_dat$d == 31] <- 31
    
    hyp_dat$Demand[hyp_dat$d == 5.5] <- 5.5
    hyp_dat$Demand[hyp_dat$d == 6.5] <- 6.5
    hyp_dat$Demand[hyp_dat$d == 7.5] <- 7.5
    hyp_dat$Demand[hyp_dat$d == 8.5] <- 8.5
    hyp_dat$Demand[hyp_dat$d == 9.5] <- 9.5
    hyp_dat$Demand[hyp_dat$d == 10.5] <- 10.5
    hyp_dat$Demand[hyp_dat$d == 11.5] <- 11.5
    hyp_dat$Demand[hyp_dat$d == 12.5] <- 12.5
    hyp_dat$Demand[hyp_dat$d == 13.5] <- 13.5
    hyp_dat$Demand[hyp_dat$d == 14.5] <- 14.5
    
    hyp_dat$`Sanctions Cost` <- NA
    
    hyp_dat$`Sanctions Cost`[hyp_dat$sc == 1 ] <- 1 
    hyp_dat$`Sanctions Cost`[hyp_dat$sc == 2 ] <- 2 
    hyp_dat$`Sanctions Cost`[hyp_dat$sc == 3 ] <- 3 
    hyp_dat$`Sanctions Cost`[hyp_dat$sc == 4 ] <- 4 
    hyp_dat$`Sanctions Cost`[hyp_dat$sc == 5 ] <- 5 
    hyp_dat$`Sanctions Cost`[hyp_dat$sc == 6 ] <- 6 
    hyp_dat$`Sanctions Cost`[hyp_dat$sc == 7 ] <- 7 
    hyp_dat$`Sanctions Cost`[hyp_dat$sc == 8 ] <- 8 
    hyp_dat$`Sanctions Cost`[hyp_dat$sc == 9 ] <- 9 
    hyp_dat$`Sanctions Cost`[hyp_dat$sc == 10] <- 10
    hyp_dat$`Sanctions Cost`[hyp_dat$sc == 11] <- 11
    hyp_dat$`Sanctions Cost`[hyp_dat$sc == 12] <- 12
    hyp_dat$`Sanctions Cost`[hyp_dat$sc == 13] <- 13
    hyp_dat$`Sanctions Cost`[hyp_dat$sc == 14] <- 14
    hyp_dat$`Sanctions Cost`[hyp_dat$sc == 15] <- 15
    hyp_dat$`Sanctions Cost`[hyp_dat$sc == 16] <- 16
    hyp_dat$`Sanctions Cost`[hyp_dat$sc == 17] <- 17
    hyp_dat$`Sanctions Cost`[hyp_dat$sc == 18] <- 18
    hyp_dat$`Sanctions Cost`[hyp_dat$sc == 19] <- 19
    hyp_dat$`Sanctions Cost`[hyp_dat$sc == 20] <- 20
    
    hyp_dat <- hyp_dat |> rename(`Sanctions Demand` =  Demand)
    
    hyp_dat$Effect_n <- NA
    
    hyp_dat$Effect_n[hyp_dat$Effect == "Demand Elasticity"] <- "Precipitant x Cost"
    hyp_dat$Effect_n[hyp_dat$Effect == "Demand Shift"] <- "Precipitant + Cost"
    
    
    pdf(file="hypothesis.pdf", height=8, width=16)
    
    ggplot(hyp_dat, aes(x = `Sanctions Demand`, y = `Sanctions Cost`, lty = V4)) +
      theme_bw() +
      geom_line(lwd = 1.5) + 
      facet_wrap(~ Effect + Effect_n) +
      theme(axis.title = element_text(size = 20)) +
      theme(axis.text = element_text(size = 20)) +
      theme(strip.text = element_text(size = 18)) +
      theme(legend.position = "none") +
      theme(axis.text.x = element_blank(),
            axis.text.y = element_blank(),
            axis.ticks = element_blank())
    
    dev.off()

# Figure 3: Precipitant, Cost, and Support for Economic Sanctions ---------

  # Baseline
    mean(spbbdat$sanctionynn)
    
  # Bolivia-Brazil (Support)
    bb1 <- summary(glm(sanctionynn ~ territory + cost, 
               family = "binomial", data = bb))
    bb2 <- summary(glm(sanctionynn ~ territory + cost + territory:cost, 
                       family = "binomial", data = bb))
    bb3 <- summary(glm(sanctionynn ~ territory + cost +
                        as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       family = "binomial", data = bb))
    bb4 <- summary(glm(sanctionynn ~ territory + cost + territory:cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       family = "binomial", data = bb))
    
    bb1mod <- glm(sanctionynn ~ territory + cost, 
                       family = "binomial", data = bb)
    bb2mod <- glm(sanctionynn ~ territory + cost + territory:cost, 
                       family = "binomial", data = bb)
    bb3mod <- glm(sanctionynn ~ territory + cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       family = "binomial", data = bb)
    bb4mod <- glm(sanctionynn ~ territory + cost + territory:cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       family = "binomial", data = bb)
    
    texreg::texreg(list(bb1mod,bb2mod,bb3mod,bb4mod),
                   stars = c(.01,.05,.1),
                   digits = 3,
                   custom.coef.map = list("territory" = "Precipitant",
                                          "cost" = "Cost",
                                          "as.factor(partisan3)1" = "Independent",
                                          "as.factor(partisan3)2" = "Republican",
                                          "k40" = "Income",
                                          "millenial" = "Millenial",
                                          "white" = "White",
                                          "college" = "College Educated",
                                          "hinf" = "High Information",
                                          "salary" = "Employed for Salary",
                                          "territory:cost" = "Precipitant x Cost",
                                          "(Intercept)" = "Constant"))  
  
  # Bolivia-Brazil (Handling)
    bb5 <- summary(lm(handlingn ~ territory + cost, data = bb))
    bb6 <- summary(lm(handlingn ~ territory + cost + territory:cost, data = bb))
    bb7 <- summary(lm(handlingn ~ territory + cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                      data = bb))
    bb8 <- summary(lm(handlingn ~ territory + cost + territory:cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                      data = bb))
    
    bb5mod <- lm(handlingn ~ territory + cost, data = bb)
    bb6mod <- lm(handlingn ~ territory + cost + territory:cost, data = bb)
    bb7mod <- lm(handlingn ~ territory + cost +
                        as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                      data = bb)
    bb8mod <- lm(handlingn ~ territory + cost + territory:cost +
                        as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                      data = bb)
    
    texreg::texreg(list(bb5mod,bb6mod,bb7mod,bb8mod),
                   stars = c(.01,.05,.1),
                   digits = 3,
                   custom.coef.map = list("territory" = "Precipitant",
                                          "cost" = "Cost",
                                          "as.factor(partisan3)1" = "Independent",
                                          "as.factor(partisan3)2" = "Republican",
                                          "k40" = "Income",
                                          "millenial" = "Millenial",
                                          "white" = "White",
                                          "college" = "College Educated",
                                          "hinf" = "High Information",
                                          "salary" = "Employed for Salary",
                                          "territory:cost" = "Precipitant x Cost",
                                          "(Intercept)" = "Constant"))  
    
  # Bolivia-Brazil (Approval)  
    bb9  <- summary(lm(approvaln ~ territory + cost, data = bb))
    bb10 <- summary(lm(approvaln ~ territory + cost + territory:cost, data = bb))
    bb11 <- summary(lm(approvaln ~ territory + cost +
                        as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                      data = bb))
    bb12 <- summary(lm(approvaln ~ territory + cost + territory:cost +
                        as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                      data = bb))
    
    bb9mod  <- lm(approvaln ~ territory + cost, data = bb)
    bb10mod <- lm(approvaln ~ territory + cost + territory:cost, data = bb)
    bb11mod <- lm(approvaln ~ territory + cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       data = bb)
    bb12mod <- lm(approvaln ~ territory + cost + territory:cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       data = bb)
    
    texreg::texreg(list(bb9mod,bb10mod,bb11mod,bb12mod),
                   stars = c(.01,.05,.1),
                   digits = 3,
                   custom.coef.map = list("territory" = "Precipitant",
                                          "cost" = "Cost",
                                          "as.factor(partisan3)1" = "Independent",
                                          "as.factor(partisan3)2" = "Republican",
                                          "k40" = "Income",
                                          "millenial" = "Millenial",
                                          "white" = "White",
                                          "college" = "College Educated",
                                          "hinf" = "High Information",
                                          "salary" = "Employed for Salary",
                                          "territory:cost" = "Precipitant x Cost",
                                          "(Intercept)" = "Constant"))  
    
  # Bolivia-Brazil (Competence)
    bb13 <- summary(lm(competencen ~ territory + cost, data = bb))
    bb14 <- summary(lm(competencen ~ territory + cost + territory:cost, data = bb))
    bb15 <- summary(lm(competencen ~ territory + cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       data = bb))
    bb16 <- summary(lm(competencen ~ territory + cost + territory:cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       data = bb))
    
    bb13mod <- summary(lm(competencen ~ territory + cost, data = bb))
    bb14mod <- summary(lm(competencen ~ territory + cost + territory:cost, data = bb))
    bb15mod <- summary(lm(competencen ~ territory + cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       data = bb))
    bb16mod <- summary(lm(competencen ~ territory + cost + territory:cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       data = bb))
    
    texreg::texreg(list(bb13mod,bb14mod,bb15mod,bb16mod),
                   stars = c(.01,.05,.1),
                   digits = 3,
                   custom.coef.map = list("territory" = "Precipitant",
                                          "cost" = "Cost",
                                          "as.factor(partisan3)1" = "Independent",
                                          "as.factor(partisan3)2" = "Republican",
                                          "k40" = "Income",
                                          "millenial" = "Millenial",
                                          "white" = "White",
                                          "college" = "College Educated",
                                          "hinf" = "High Information",
                                          "salary" = "Employed for Salary",
                                          "territory:cost" = "Precipitant x Cost",
                                          "(Intercept)" = "Constant")) 
    
    
  # Spain-Portugal (Support)
    sp1 <- summary(glm(sanctionynn ~ territory + cost, 
                       family = "binomial", data = sp))
    sp2 <- summary(glm(sanctionynn ~ territory + cost + territory:cost, 
                       family = "binomial", data = sp))
    sp3 <- summary(glm(sanctionynn ~ territory + cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       family = "binomial", data = sp))
    sp4 <- summary(glm(sanctionynn ~ territory + cost + territory:cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       family = "binomial", data = sp))
    
    sp1mod <- glm(sanctionynn ~ territory + cost, 
                       family = "binomial", data = sp)
    sp2mod <- glm(sanctionynn ~ territory + cost + territory:cost, 
                       family = "binomial", data = sp)
    sp3mod <- glm(sanctionynn ~ territory + cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       family = "binomial", data = sp)
    sp4mod <- glm(sanctionynn ~ territory + cost + territory:cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       family = "binomial", data = sp)
    
    
    texreg::texreg(list(sp1mod,sp2mod,sp3mod,sp4mod),
                   stars = c(.01,.05,.1),
                   digits = 3,
                   custom.coef.map = list("territory" = "Precipitant",
                                          "cost" = "Cost",
                                          "as.factor(partisan3)1" = "Independent",
                                          "as.factor(partisan3)2" = "Republican",
                                          "k40" = "Income",
                                          "millenial" = "Millenial",
                                          "white" = "White",
                                          "college" = "College Educated",
                                          "hinf" = "High Information",
                                          "salary" = "Employed for Salary",
                                          "territory:cost" = "Precipitant x Cost",
                                          "(Intercept)" = "Constant"))  
    
  # Spain-Portugal (Handling)
    sp5 <- summary(lm(handlingn ~ territory + cost, data = sp))
    sp6 <- summary(lm(handlingn ~ territory + cost + territory:cost, data = sp))
    sp7 <- summary(lm(handlingn ~ territory + cost +
                        as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                      data = sp))
    sp8 <- summary(lm(handlingn ~ territory + cost + territory:cost +
                        as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                      data = sp))
    
    sp5mod <- lm(handlingn ~ territory + cost, data = sp)
    sp6mod <- lm(handlingn ~ territory + cost + territory:cost, data = sp)
    sp7mod <- lm(handlingn ~ territory + cost +
                        as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                      data = sp)
    sp8mod <- lm(handlingn ~ territory + cost + territory:cost +
                        as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                      data = sp)
    
    texreg::texreg(list(sp5mod,sp6mod,sp7mod,sp8mod),
                   stars = c(.01,.05,.1),
                   digits = 3,
                   custom.coef.map = list("territory" = "Precipitant",
                                          "cost" = "Cost",
                                          "as.factor(partisan3)1" = "Independent",
                                          "as.factor(partisan3)2" = "Republican",
                                          "k40" = "Income",
                                          "millenial" = "Millenial",
                                          "white" = "White",
                                          "college" = "College Educated",
                                          "hinf" = "High Information",
                                          "salary" = "Employed for Salary",
                                          "territory:cost" = "Precipitant x Cost",
                                          "(Intercept)" = "Constant"))  
    
  # Spain-Portugal (Approval)  
    sp9  <- summary(lm(approvaln ~ territory + cost, data = sp))
    sp10 <- summary(lm(approvaln ~ territory + cost + territory:cost, data = sp))
    sp11 <- summary(lm(approvaln ~ territory + cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       data = sp))
    sp12 <- summary(lm(approvaln ~ territory + cost + territory:cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       data = sp))
    
    sp9mod  <-lm(approvaln ~ territory + cost, data = sp)
    sp10mod <-lm(approvaln ~ territory + cost + territory:cost, data = sp)
    sp11mod <-lm(approvaln ~ territory + cost +
                as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
              data = sp)
    sp12mod <-lm(approvaln ~ territory + cost + territory:cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       data = sp)
    
    texreg::texreg(list(sp9mod,sp10mod,sp11mod,sp12mod),
                   stars = c(.01,.05,.1),
                   digits = 3,
                   custom.coef.map = list("territory" = "Precipitant",
                                          "cost" = "Cost",
                                          "as.factor(partisan3)1" = "Independent",
                                          "as.factor(partisan3)2" = "Republican",
                                          "k40" = "Income",
                                          "millenial" = "Millenial",
                                          "white" = "White",
                                          "college" = "College Educated",
                                          "hinf" = "High Information",
                                          "salary" = "Employed for Salary",
                                          "territory:cost" = "Precipitant x Cost",
                                          "(Intercept)" = "Constant"))  
    
  # Spain-Portugal (Competence)  
    sp13 <- summary(lm(competencen ~ territory + cost, data = sp))
    sp14 <- summary(lm(competencen ~ territory + cost + territory:cost, data = sp))
    sp15 <- summary(lm(competencen ~ territory + cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       data = sp))
    sp16 <- summary(lm(competencen ~ territory + cost + territory:cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       data = sp))
    
    sp13mod <- lm(competencen ~ territory + cost, data = sp)
    sp14mod <- lm(competencen ~ territory + cost + territory:cost, data = sp)
    sp15mod <- lm(competencen ~ territory + cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       data = sp)
    sp16mod <- lm(competencen ~ territory + cost + territory:cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       data = sp)
    
    texreg::texreg(list(sp13mod,sp14mod,sp15mod,sp16mod),
                   stars = c(.01,.05,.1),
                   digits = 3,
                   custom.coef.map = list("territory" = "Precipitant",
                                          "cost" = "Cost",
                                          "as.factor(partisan3)1" = "Independent",
                                          "as.factor(partisan3)2" = "Republican",
                                          "k40" = "Income",
                                          "millenial" = "Millenial",
                                          "white" = "White",
                                          "college" = "College Educated",
                                          "hinf" = "High Information",
                                          "salary" = "Employed for Salary",
                                          "territory:cost" = "Precipitant x Cost",
                                          "(Intercept)" = "Constant"))  
    
  # Pooled (Support)
    pool1 <- summary(glm(sanctionynn ~ territory + cost, 
                         family = "binomial", data = spbbdat))
    pool2 <- summary(glm(sanctionynn ~ territory + cost + territory:cost, 
                         family = "binomial", data = spbbdat))
    pool3 <- summary(glm(sanctionynn ~ territory + cost +
                           as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                         family = "binomial", data = spbbdat))
    pool4 <- summary(glm(sanctionynn ~ territory + cost + territory:cost +
                           as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                         family = "binomial", data = spbbdat))
    
    pool1mod <- glm(sanctionynn ~ territory + cost, 
                 family = "binomial", data = spbbdat)
    pool2mod <- glm(sanctionynn ~ territory + cost + territory:cost, 
                 family = "binomial", data = spbbdat)
    pool3mod <- glm(sanctionynn ~ territory + cost +
                   as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                 family = "binomial", data = spbbdat)
    pool4mod <- glm(sanctionynn ~ territory + cost + territory:cost +
                           as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                         family = "binomial", data = spbbdat)
    
    
    texreg::texreg(list(pool1mod,pool2mod,pool3mod,pool4mod),
                   stars = c(.01,.05,.1),
                   digits = 3,
                   custom.coef.map = list("territory" = "Precipitant",
                                          "cost" = "Cost",
                                          "as.factor(partisan3)1" = "Independent",
                                          "as.factor(partisan3)2" = "Republican",
                                          "k40" = "Income",
                                          "millenial" = "Millenial",
                                          "white" = "White",
                                          "college" = "College Educated",
                                          "hinf" = "High Information",
                                          "salary" = "Employed for Salary",
                                          "territory:cost" = "Precipitant x Cost",
                                          "(Intercept)" = "Constant")) 
    
  # Pooled (Handling)  
    pool5 <- summary(lm(handlingn ~ territory + cost, data = spbbdat))
    pool6 <- summary(lm(handlingn ~ territory + cost + territory:cost, data = spbbdat))
    pool7 <- summary(lm(handlingn ~ territory + cost +
                          as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                        data = spbbdat))
    pool8 <- summary(lm(handlingn ~ territory + cost + territory:cost +
                          as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                        data = spbbdat))
    
    
    pool5mod <- lm(handlingn ~ territory + cost, data = spbbdat)
    pool6mod <- lm(handlingn ~ territory + cost + territory:cost, data = spbbdat)
    pool7mod <- lm(handlingn ~ territory + cost +
                  as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                data = spbbdat)
    pool8mod <- lm(handlingn ~ territory + cost + territory:cost +
                          as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                        data = spbbdat)
    
    texreg::texreg(list(pool5mod,pool6mod,pool7mod,pool8mod),
                   stars = c(.01,.05,.1),
                   digits = 3,
                   custom.coef.map = list("territory" = "Precipitant",
                                          "cost" = "Cost",
                                          "as.factor(partisan3)1" = "Independent",
                                          "as.factor(partisan3)2" = "Republican",
                                          "k40" = "Income",
                                          "millenial" = "Millenial",
                                          "white" = "White",
                                          "college" = "College Educated",
                                          "hinf" = "High Information",
                                          "salary" = "Employed for Salary",
                                          "territory:cost" = "Precipitant x Cost",
                                          "(Intercept)" = "Constant")) 
    
  # Pooled (Presidential Approval)  
    pool9  <- summary(lm(approvaln ~ territory + cost, data = spbbdat))
    pool10 <- summary(lm(approvaln ~ territory + cost + territory:cost, data = spbbdat))
    pool11 <- summary(lm(approvaln ~ territory + cost +
                           as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                         data = spbbdat))
    pool12 <- summary(lm(approvaln ~ territory + cost + territory:cost +
                           as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                         data = spbbdat))
    
    
    pool9mod  <- summary(lm(approvaln ~ territory + cost, data = spbbdat))
    pool10mod <- summary(lm(approvaln ~ territory + cost + territory:cost, data = spbbdat))
    pool11mod <- summary(lm(approvaln ~ territory + cost +
                           as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                         data = spbbdat))
    pool12mod <- summary(lm(approvaln ~ territory + cost + territory:cost +
                           as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                         data = spbbdat))
    
    
    texreg::texreg(list(pool9mod,pool10mod,pool11mod,pool12mod),
                   stars = c(.01,.05,.1),
                   digits = 3,
                   custom.coef.map = list("territory" = "Precipitant",
                                          "cost" = "Cost",
                                          "as.factor(partisan3)1" = "Independent",
                                          "as.factor(partisan3)2" = "Republican",
                                          "k40" = "Income",
                                          "millenial" = "Millenial",
                                          "white" = "White",
                                          "college" = "College Educated",
                                          "hinf" = "High Information",
                                          "salary" = "Employed for Salary",
                                          "territory:cost" = "Precipitant x Cost",
                                          "(Intercept)" = "Constant")) 
    
  # Pooled (Competence)  
    pool13 <- summary(lm(competencen ~ territory + cost, data = spbbdat))
    pool14 <- summary(lm(competencen ~ territory + cost + territory:cost, data = spbbdat))
    pool15 <- summary(lm(competencen ~ territory + cost +
                           as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                         data = spbbdat))
    pool16 <- summary(lm(competencen ~ territory + cost + territory:cost +
                           as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                         data = spbbdat))
    
    pool13mod <- lm(competencen ~ territory + cost, data = spbbdat)
    pool14mod <- lm(competencen ~ territory + cost + territory:cost, data = spbbdat)
    pool15mod <- lm(competencen ~ territory + cost +
                   as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                 data = spbbdat)
    pool16mod <- lm(competencen ~ territory + cost + territory:cost +
                           as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                         data = spbbdat)
    
    texreg::texreg(list(pool13mod,pool14mod,pool15mod,pool16mod),
                   stars = c(.01,.05,.1),
                   digits = 3,
                   custom.coef.map = list("territory" = "Precipitant",
                                          "cost" = "Cost",
                                          "as.factor(partisan3)1" = "Independent",
                                          "as.factor(partisan3)2" = "Republican",
                                          "k40" = "Income",
                                          "millenial" = "Millenial",
                                          "white" = "White",
                                          "college" = "College Educated",
                                          "hinf" = "High Information",
                                          "salary" = "Employed for Salary",
                                          "territory:cost" = "Precipitant x Cost",
                                          "(Intercept)" = "Constant")) 
    
  # Competence as a function of preference
    summary(lm(competencen ~ sanctionynn + territory + cost + territory:cost +
               as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
               data = spbbdat))
    
  # Set Up Results Matrix  
    Countries <- c("Spain-Portugal","Brazil-Bolivia","Pooled")
    Effect <- c("Conflict","Cost","Conflict x Cost")
    
    results_sanc <- as_tibble(expand.grid("Countries" = Countries,
                              "Effect" = Effect))
    results_sanc
    
  # Harvest Effects
    results_sanc$`Odds Ratios` <- NA

    results_sanc$`Odds Ratios`[results_sanc$Countries == "Spain-Portugal" & results_sanc$Effect == "Conflict"] <- round((exp(sp3$coefficients["territory","Estimate"]) - 1)*100,1)
    results_sanc$`Odds Ratios`[results_sanc$Countries == "Spain-Portugal" & results_sanc$Effect == "Cost"]     <- round((exp(sp3$coefficients["cost","Estimate"]) - 1)*100,1)
    
    results_sanc$`Odds Ratios`[results_sanc$Countries == "Brazil-Bolivia" & results_sanc$Effect == "Conflict"] <- round((exp(bb3$coefficients["territory","Estimate"]) - 1)*100,1)
    results_sanc$`Odds Ratios`[results_sanc$Countries == "Brazil-Bolivia" & results_sanc$Effect == "Cost"]     <- round((exp(bb3$coefficients["cost","Estimate"]) - 1)*100,1)
          
    results_sanc$`Odds Ratios`[results_sanc$Countries == "Pooled" & results_sanc$Effect == "Conflict"] <- round((exp(pool3$coefficients["territory","Estimate"]) - 1)*100,1)
    results_sanc$`Odds Ratios`[results_sanc$Countries == "Pooled" & results_sanc$Effect == "Cost"]     <- round((exp(pool3$coefficients["cost","Estimate"]) - 1)*100,1)
    
    results_sanc$`Odds Ratios`[results_sanc$Countries == "Spain-Portugal" & results_sanc$Effect == "Conflict x Cost"] <- round((exp(sp4$coefficients["territory:cost","Estimate"]) - 1)*100,1)
    results_sanc$`Odds Ratios`[results_sanc$Countries == "Brazil-Bolivia" & results_sanc$Effect == "Conflict x Cost"] <- round((exp(bb4$coefficients["territory:cost","Estimate"]) - 1)*100,1)
    results_sanc$`Odds Ratios`[results_sanc$Countries == "Pooled" & results_sanc$Effect == "Conflict x Cost"] <- round((exp(pool4$coefficients["territory:cost","Estimate"]) - 1)*100,1)
    
    results_sanc$se <- NA
    results_sanc$se[results_sanc$Countries == "Spain-Portugal" & results_sanc$Effect == "Conflict"] <- round((exp(sp3$coefficients["territory","Std. Error"]) - 1)*100,1)
    results_sanc$se[results_sanc$Countries == "Spain-Portugal" & results_sanc$Effect == "Cost"]     <- round((exp(sp3$coefficients["cost","Std. Error"]) - 1)*100,1)
    results_sanc$se[results_sanc$Countries == "Brazil-Bolivia" & results_sanc$Effect == "Conflict"] <- round((exp(bb3$coefficients["territory","Std. Error"]) - 1)*100,1)
    results_sanc$se[results_sanc$Countries == "Brazil-Bolivia" & results_sanc$Effect == "Cost"]     <- round((exp(bb3$coefficients["cost","Std. Error"]) - 1)*100,1)
    
    results_sanc$se[results_sanc$Countries == "Pooled" & results_sanc$Effect == "Conflict"] <- round((exp(pool3$coefficients["territory","Std. Error"]) - 1)*100,1)
    results_sanc$se[results_sanc$Countries == "Pooled" & results_sanc$Effect == "Cost"]     <- round((exp(pool3$coefficients["cost","Std. Error"]) - 1)*100,1)
    results_sanc$se[results_sanc$Countries == "Spain-Portugal" & results_sanc$Effect == "Conflict x Cost"] <- round((exp(sp4$coefficients["territory:cost","Std. Error"]) - 1)*100,1)
    results_sanc$se[results_sanc$Countries == "Brazil-Bolivia" & results_sanc$Effect == "Conflict x Cost"] <- round((exp(bb4$coefficients["territory:cost","Std. Error"]) - 1)*100,1)
    results_sanc$se[results_sanc$Countries == "Pooled" & results_sanc$Effect == "Conflict x Cost"] <- round((exp(pool4$coefficients["territory:cost","Std. Error"]) - 1)*100,1)
    
    results_sanc$lb <- NA
    
    results_sanc$lb <- results_sanc$`Odds Ratios` - results_sanc$se*1.68
    results_sanc$ub <- results_sanc$`Odds Ratios` + results_sanc$se*1.68
    
    results_sanc$Effect_n <- NA
    results_sanc$Effect_n[results_sanc$Effect == "Conflict"] <- "C"
    results_sanc$Effect_n[results_sanc$Effect == "Cost"] <- "B"
    results_sanc$Effect_n[results_sanc$Effect == "Conflict x Cost"] <- "A"
    
    results_sanc
    
    #setwd("../Plots/")
    
    pdf(file="effectsrr.pdf", height=8, width=20)
    
      ggplot(results_sanc, aes(x = `Odds Ratios`, y = Effect_n, color = Countries)) + 
        theme_bw() +
        geom_point(size = 5, position = position_dodge2(w = 0.75), shape = 15) +
        geom_linerange(aes(xmin = lb, xmax = ub), 
                       position = position_dodge2(w = 0.75),
                       lwd = 1.5) +
        scale_color_manual(values=c("lightgray", "darkgray", "black")) +
        labs(x = "\n % Change in Baseline Support for Sanctions") +
        scale_y_discrete(name = "Hypothesized Effects \n",
                         labels = c("C" = "Precipitant =\n Invasion", "B" = "Costly \n Sanctions", "A" = "Precipitant = Invasion \n x Costly Sanctions")) +
        theme(axis.title = element_text(size = 24)) +
        theme(axis.text = element_text(size = 20)) +
        theme(legend.title = element_text(size = 20, face = "bold"),
              legend.text = element_text(size = 20),
              legend.title.align=0.5) +
        theme(axis.text.y = element_text(angle=90, hjust = .5)) +
        geom_vline(xintercept = 0, lwd = 1.2, lty = 2) +
        xlim(-275,275)
    
    dev.off()
    
  # Marginal Effects

  # Partisan Variables
    spbbdat$dem <- NA
    
    spbbdat$dem[spbbdat$partisan3 == 0] <- 1
    spbbdat$dem[spbbdat$partisan3 == 1] <- 0
    spbbdat$dem[spbbdat$partisan3 == 2] <- 0
    
    spbbdat$ind <- NA
    
    spbbdat$ind[spbbdat$partisan3 == 0] <- 0
    spbbdat$ind[spbbdat$partisan3 == 1] <- 1
    spbbdat$ind[spbbdat$partisan3 == 2] <- 0
    
    spbbdat$gop <- NA
    
    spbbdat$gop[spbbdat$partisan3 == 0] <- 0
    spbbdat$gop[spbbdat$partisan3 == 1] <- 0
    spbbdat$gop[spbbdat$partisan3 == 2] <- 1
  
    table(spbbdat$partisan3)
    table(spbbdat$dem)
    table(spbbdat$ind)
    table(spbbdat$gop)
      
  # Modify data
    spbbdat$cost_f <- as.factor(spbbdat$cost)
    spbbdat$territory_f <- as.factor(spbbdat$territory)
    
    summary(pool4mod <- glm(sanctionynn ~ territory + cost + territory:cost +
                    as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                    family = "binomial", data = spbbdat))
    
    summary(sp4mod <- glm(sanctionynn ~ territory + cost + territory:cost +
                    as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                  family = "binomial", data = sp))
    
    summary(bb4mod <- glm(sanctionynn ~ territory + cost + territory:cost +
                    as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                  family = "binomial", data = bb))
    
    
   mfx_cost_pooled <-  interplot(m = pool4mod, var1 = "cost", var2 = "territory", ci = .9) +
      labs(
        title = "Marginal Effects of Cost at Different Levels of Territory",
        x = "Territory Levels",
        y = "Marginal Effect of Cost"
      ) 
   
   mfx_cost_pooled <- as_tibble(mfx_cost_pooled$data)
   
   mfx_cost_pooled
   
   mfx_cost_pooled$Moderator <- NA
   mfx_cost_pooled$Moderator[mfx_cost_pooled$fake == 0] <- "Trade"
   mfx_cost_pooled$Moderator[mfx_cost_pooled$fake == 1] <- "Invasion"
    
   mfx_prec_pooled <- interplot(m = pool4mod, var1 = "territory", var2 = "cost", ci = .9) +
      labs(
        title = "Marginal Effects of Cost at Different Levels of Territory",
        x = "Territory Levels",
        y = "Marginal Effect of Cost"
      ) 
   
   mfx_prec_pooled <- as_tibble(mfx_prec_pooled$data)
   
   mfx_prec_pooled$Moderator <- NA
   mfx_prec_pooled$Moderator[mfx_prec_pooled$fake == 0] <- "Not Costly"
   mfx_prec_pooled$Moderator[mfx_prec_pooled$fake == 1] <- "Costly"
   
   mfx_prec_pooled
   
   mfx_cost_sp <-  interplot(m = sp4mod, var1 = "cost", var2 = "territory", ci = .9) +
     labs(
       title = "Marginal Effects of Cost at Different Levels of Territory",
       x = "Territory Levels",
       y = "Marginal Effect of Cost"
     ) 
   
   mfx_cost_sp <- as_tibble(mfx_cost_sp$data)
   
   mfx_cost_sp
   
   mfx_cost_sp$Moderator <- NA
   mfx_cost_sp$Moderator[mfx_cost_sp$fake == 0] <- "Trade"
   mfx_cost_sp$Moderator[mfx_cost_sp$fake == 1] <- "Invasion"
   
   mfx_prec_sp <- interplot(m = sp4mod, var1 = "territory", var2 = "cost", ci = .9) +
     labs(
       title = "Marginal Effects of Cost at Different Levels of Territory",
       x = "Territory Levels",
       y = "Marginal Effect of Cost"
     ) 
   
   mfx_prec_sp <- as_tibble(mfx_prec_sp$data)
   
   mfx_prec_sp$Moderator <- NA
   mfx_prec_sp$Moderator[mfx_prec_sp$fake == 0] <- "Not Costly"
   mfx_prec_sp$Moderator[mfx_prec_sp$fake == 1] <- "Costly"
   
   mfx_prec_sp
   
   mfx_cost_bb <-  interplot(m = bb4mod, var1 = "cost", var2 = "territory", ci = .9) +
     labs(
       title = "Marginal Effects of Cost at Different Levels of Territory",
       x = "Territory Levels",
       y = "Marginal Effect of Cost"
     ) 
   
   mfx_cost_bb <- as_tibble(mfx_cost_bb$data)
   
   mfx_cost_bb
   
   mfx_cost_bb$Moderator <- NA
   mfx_cost_bb$Moderator[mfx_cost_bb$fake == 0] <- "Trade"
   mfx_cost_bb$Moderator[mfx_cost_bb$fake == 1] <- "Invasion"
   
   mfx_cost_pooled$Sample <- "Pooled"
   mfx_cost_sp$Sample <- "Spain-Portugal"
   mfx_cost_bb$Sample <- "Brazil-Bolivia"
   
   cost_mfx <- rbind(mfx_cost_bb,
                     mfx_cost_sp,
                     mfx_cost_pooled)
   
   mfx_prec_bb <- interplot(m = bb4mod, var1 = "territory", var2 = "cost", ci = .9) +
     labs(
       title = "Marginal Effects of Cost at Different Levels of Territory",
       x = "Territory Levels",
       y = "Marginal Effect of Cost"
     ) 
   
   mfx_prec_bb <- as_tibble(mfx_prec_bb$data)
   
   mfx_prec_bb$Moderator <- NA
   mfx_prec_bb$Moderator[mfx_prec_bb$fake == 0] <- "Not Costly"
   mfx_prec_bb$Moderator[mfx_prec_bb$fake == 1] <- "Costly"
   
   mfx_prec_bb
   
   mfx_prec_pooled$Sample <- "Pooled"
   mfx_prec_sp$Sample <- "Spain-Portugal"
   mfx_prec_bb$Sample <- "Brazil-Bolivia"
   
   prec_mfx <- rbind(mfx_prec_pooled,
                     mfx_prec_sp,
                     mfx_prec_bb)
   
   prec_mfx
   
   pdf(file="mfxprec.pdf", height=8, width=10)
   
   ggplot(prec_mfx, aes(x = Moderator, y = coef1, color = Sample)) +
     geom_point(size = 5, shape = 15, position = position_dodge(width = 0.5)) +
     #geom_errorbar(aes(ymin = lb, ymax = ub), width = 0.25, 
     #position = position_dodge(width = 0.5)) +
     geom_linerange(aes(ymin = lb, ymax = ub), 
                    position = position_dodge(width = 0.5),
                    lwd = 1.5) +
     theme_bw() +
     ylim(-2,2) +
     geom_hline(yintercept = 0, lty = 2, lwd = 2) +
     scale_color_manual(values = c("darkgray", "black", "lightgray")) +
     labs(x = "Cost of Sanctions", y = "Effect of Precipitant", color = "Sample") +
     theme(axis.title = element_text(size = 20)) +
     theme(axis.text = element_text(size = 20)) +
     theme(legend.title = element_text(size = 20, face = "bold"),
           legend.text = element_text(size = 20),
           strip.text = element_text(size = 20),
           legend.title.align=0.5,
           legend.position = "none") +
     theme(axis.text=element_text(size=18),
           axis.title=element_text(size=24),
           plot.margin = margin(0.5,0.5,1.1,1.1,"cm"),
           axis.title.x = element_text(vjust = -2),
           axis.title.y = element_text(vjust = 4))
   
   dev.off()
   
   pdf(file="mfxcost.pdf", height=8, width=10)
   
   ggplot(cost_mfx, aes(x = Moderator, y = coef1, color = Sample)) +
     geom_point(size = 5, shape = 15, position = position_dodge(width = 0.5)) +
     #geom_errorbar(aes(ymin = lb, ymax = ub), width = 0.25, 
     #               position = position_dodge(width = 0.5)) +
     geom_linerange(aes(ymin = lb, ymax = ub), 
                    position = position_dodge(width = 0.5),
                    lwd = 1.5) +
     theme_bw() +
     ylim(-2,2) +
     geom_hline(yintercept = 0, lty = 2, lwd = 2) +
     scale_color_manual(values = c("darkgray", "black", "lightgray")) +
     labs(x = "Precipitant", y = "Effect of Cost", color = "Sample") +
     theme(axis.title = element_text(size = 20)) +
     theme(axis.text = element_text(size = 20)) +
     theme(legend.title = element_text(size = 20, face = "bold"),
           legend.text = element_text(size = 20),
           strip.text = element_text(size = 20),
           legend.title.align=0.5,
           legend.position = "none") +
     theme(axis.text=element_text(size=18),
           axis.title=element_text(size=24),
           plot.margin = margin(0.5,0.5,1.1,1.1,"cm"),
           axis.title.x = element_text(vjust = -2),
           axis.title.y = element_text(vjust = 4))
   
   dev.off()
     
# Figure 4: The Political Consequences of Precipitant Severity and --------

  # Set Up Results  
    Countries <- c("Spain-Portugal","Brazil-Bolivia","Pooled")
    Effect <- c("Conflict","Cost","Conflict x Cost")
    
    results_sanc_pols <- as_tibble(expand.grid("Countries" = Countries,
                                          "Effect" = Effect))
    results_sanc_pols
    
    
  # Handling
    results_sanc_pols$Handling <- NA
    
    results_sanc_pols$Handling[results_sanc_pols$Countries == "Spain-Portugal" & results_sanc_pols$Effect == "Conflict"] <- round(sp7$coefficients["territory","Estimate"],2)
    results_sanc_pols$Handling[results_sanc_pols$Countries == "Spain-Portugal" & results_sanc_pols$Effect == "Cost"]     <- round(sp7$coefficients["cost","Estimate"],2)
    
    results_sanc_pols$Handling[results_sanc_pols$Countries == "Brazil-Bolivia" & results_sanc_pols$Effect == "Conflict"] <- round(bb7$coefficients["territory","Estimate"],2)
    results_sanc_pols$Handling[results_sanc_pols$Countries == "Brazil-Bolivia" & results_sanc_pols$Effect == "Cost"]     <- round(bb7$coefficients["cost","Estimate"],2)
    
    results_sanc_pols$Handling[results_sanc_pols$Countries == "Pooled" & results_sanc_pols$Effect == "Conflict"] <- round(pool7$coefficients["territory","Estimate"],2)
    results_sanc_pols$Handling[results_sanc_pols$Countries == "Pooled" & results_sanc_pols$Effect == "Cost"]     <- round(pool7$coefficients["cost","Estimate"],2)
    
    results_sanc_pols$Handling[results_sanc_pols$Countries == "Spain-Portugal" & results_sanc_pols$Effect == "Conflict x Cost"] <- round(sp8$coefficients["territory:cost","Estimate"],2)
    results_sanc_pols$Handling[results_sanc_pols$Countries == "Brazil-Bolivia" & results_sanc_pols$Effect == "Conflict x Cost"] <- round(bb8$coefficients["territory:cost","Estimate"],2)
    results_sanc_pols$Handling[results_sanc_pols$Countries == "Pooled" & results_sanc_pols$Effect == "Conflict x Cost"] <- round(pool8$coefficients["territory:cost","Estimate"],2)
    
    results_sanc_pols$h_se <- NA
    results_sanc_pols$h_se[results_sanc_pols$Countries == "Spain-Portugal" & results_sanc_pols$Effect == "Conflict"] <- round(sp7$coefficients["territory","Std. Error"],2)
    results_sanc_pols$h_se[results_sanc_pols$Countries == "Spain-Portugal" & results_sanc_pols$Effect == "Cost"]     <- round(sp7$coefficients["cost","Std. Error"],2)
    results_sanc_pols$h_se[results_sanc_pols$Countries == "Brazil-Bolivia" & results_sanc_pols$Effect == "Conflict"] <- round(bb7$coefficients["territory","Std. Error"],2)
    results_sanc_pols$h_se[results_sanc_pols$Countries == "Brazil-Bolivia" & results_sanc_pols$Effect == "Cost"]     <- round(bb7$coefficients["cost","Std. Error"],2)
    
    results_sanc_pols$h_se[results_sanc_pols$Countries == "Pooled" & results_sanc_pols$Effect == "Conflict"] <- round(pool7$coefficients["territory","Std. Error"],2)
    results_sanc_pols$h_se[results_sanc_pols$Countries == "Pooled" & results_sanc_pols$Effect == "Cost"]     <- round(pool7$coefficients["cost","Std. Error"],2)
    results_sanc_pols$h_se[results_sanc_pols$Countries == "Spain-Portugal" & results_sanc_pols$Effect == "Conflict x Cost"] <- round(sp8$coefficients["territory:cost","Std. Error"],2)
    results_sanc_pols$h_se[results_sanc_pols$Countries == "Brazil-Bolivia" & results_sanc_pols$Effect == "Conflict x Cost"] <- round(bb8$coefficients["territory:cost","Std. Error"],2)
    results_sanc_pols$h_se[results_sanc_pols$Countries == "Pooled" & results_sanc_pols$Effect == "Conflict x Cost"] <- round(pool8$coefficients["territory:cost","Std. Error"],2)
    
    results_sanc_pols$lb_h <- results_sanc_pols$Handling - results_sanc_pols$h_se*1.68
    results_sanc_pols$ub_h <- results_sanc_pols$Handling + results_sanc_pols$h_se*1.68
    
    results_sanc_pols
    
    results_sanc_pols$Effect_h <- NA
    results_sanc_pols$Effect_h[results_sanc_pols$Effect == "Conflict"] <- "C"
    results_sanc_pols$Effect_h[results_sanc_pols$Effect == "Cost"] <- "B"
    results_sanc_pols$Effect_h[results_sanc_pols$Effect == "Conflict x Cost"] <- "A"
    
  # Approval
    results_sanc_pols$Approval <- NA
    
    results_sanc_pols$Approval[results_sanc_pols$Countries == "Spain-Portugal" & results_sanc_pols$Effect == "Conflict"] <- round(sp11$coefficients["territory","Estimate"],2)
    results_sanc_pols$Approval[results_sanc_pols$Countries == "Spain-Portugal" & results_sanc_pols$Effect == "Cost"]     <- round(sp11$coefficients["cost","Estimate"],2)
    
    results_sanc_pols$Approval[results_sanc_pols$Countries == "Brazil-Bolivia" & results_sanc_pols$Effect == "Conflict"] <- round(bb11$coefficients["territory","Estimate"],2)
    results_sanc_pols$Approval[results_sanc_pols$Countries == "Brazil-Bolivia" & results_sanc_pols$Effect == "Cost"]     <- round(bb11$coefficients["cost","Estimate"],2)
    
    results_sanc_pols$Approval[results_sanc_pols$Countries == "Pooled" & results_sanc_pols$Effect == "Conflict"] <- round(pool11$coefficients["territory","Estimate"],2)
    results_sanc_pols$Approval[results_sanc_pols$Countries == "Pooled" & results_sanc_pols$Effect == "Cost"]     <- round(pool11$coefficients["cost","Estimate"],2)
    
    results_sanc_pols$Approval[results_sanc_pols$Countries == "Spain-Portugal" & results_sanc_pols$Effect == "Conflict x Cost"] <- round(sp12$coefficients["territory:cost","Estimate"],2)
    results_sanc_pols$Approval[results_sanc_pols$Countries == "Brazil-Bolivia" & results_sanc_pols$Effect == "Conflict x Cost"] <- round(bb12$coefficients["territory:cost","Estimate"],2)
    results_sanc_pols$Approval[results_sanc_pols$Countries == "Pooled" & results_sanc_pols$Effect == "Conflict x Cost"] <- round(pool12$coefficients["territory:cost","Estimate"],2)
    
    results_sanc_pols$a_se <- NA
    results_sanc_pols$a_se[results_sanc_pols$Countries == "Spain-Portugal" & results_sanc_pols$Effect == "Conflict"] <- round(sp11$coefficients["territory","Std. Error"],2)
    results_sanc_pols$a_se[results_sanc_pols$Countries == "Spain-Portugal" & results_sanc_pols$Effect == "Cost"]     <- round(sp11$coefficients["cost","Std. Error"],2)
    results_sanc_pols$a_se[results_sanc_pols$Countries == "Brazil-Bolivia" & results_sanc_pols$Effect == "Conflict"] <- round(bb11$coefficients["territory","Std. Error"],2)
    results_sanc_pols$a_se[results_sanc_pols$Countries == "Brazil-Bolivia" & results_sanc_pols$Effect == "Cost"]     <- round(bb11$coefficients["cost","Std. Error"],2)
    
    results_sanc_pols$a_se[results_sanc_pols$Countries == "Pooled" & results_sanc_pols$Effect == "Conflict"] <- round(pool11$coefficients["territory","Std. Error"],2)
    results_sanc_pols$a_se[results_sanc_pols$Countries == "Pooled" & results_sanc_pols$Effect == "Cost"]     <- round(pool11$coefficients["cost","Std. Error"],2)
    results_sanc_pols$a_se[results_sanc_pols$Countries == "Spain-Portugal" & results_sanc_pols$Effect == "Conflict x Cost"] <- round(sp12$coefficients["territory:cost","Std. Error"],2)
    results_sanc_pols$a_se[results_sanc_pols$Countries == "Brazil-Bolivia" & results_sanc_pols$Effect == "Conflict x Cost"] <- round(bb12$coefficients["territory:cost","Std. Error"],2)
    results_sanc_pols$a_se[results_sanc_pols$Countries == "Pooled" & results_sanc_pols$Effect == "Conflict x Cost"] <- round(pool12$coefficients["territory:cost","Std. Error"],2)
    
    results_sanc_pols$lb_a <- results_sanc_pols$Approval - results_sanc_pols$a_se*1.612
    results_sanc_pols$ub_a <- results_sanc_pols$Approval + results_sanc_pols$a_se*1.612
    
    results_sanc_pols
    
    results_sanc_pols$Effect_a <- NA
    results_sanc_pols$Effect_a[results_sanc_pols$Effect == "Conflict"] <- "C"
    results_sanc_pols$Effect_a[results_sanc_pols$Effect == "Cost"] <- "B"
    results_sanc_pols$Effect_a[results_sanc_pols$Effect == "Conflict x Cost"] <- "A"
    
  # Competence
    results_sanc_pols$Competence <- NA
    
    results_sanc_pols$Competence[results_sanc_pols$Countries == "Spain-Portugal" & results_sanc_pols$Effect == "Conflict"] <- round(sp15$coefficients["territory","Estimate"],2)
    results_sanc_pols$Competence[results_sanc_pols$Countries == "Spain-Portugal" & results_sanc_pols$Effect == "Cost"]     <- round(sp15$coefficients["cost","Estimate"],2)
    
    results_sanc_pols$Competence[results_sanc_pols$Countries == "Brazil-Bolivia" & results_sanc_pols$Effect == "Conflict"] <- round(bb15$coefficients["territory","Estimate"],2)
    results_sanc_pols$Competence[results_sanc_pols$Countries == "Brazil-Bolivia" & results_sanc_pols$Effect == "Cost"]     <- round(bb15$coefficients["cost","Estimate"],2)
    
    results_sanc_pols$Competence[results_sanc_pols$Countries == "Pooled" & results_sanc_pols$Effect == "Conflict"] <- round(pool15$coefficients["territory","Estimate"],2)
    results_sanc_pols$Competence[results_sanc_pols$Countries == "Pooled" & results_sanc_pols$Effect == "Cost"]     <- round(pool15$coefficients["cost","Estimate"],2)
    
    results_sanc_pols$Competence[results_sanc_pols$Countries == "Spain-Portugal" & results_sanc_pols$Effect == "Conflict x Cost"] <- round(sp16$coefficients["territory:cost","Estimate"],2)
    results_sanc_pols$Competence[results_sanc_pols$Countries == "Brazil-Bolivia" & results_sanc_pols$Effect == "Conflict x Cost"] <- round(bb16$coefficients["territory:cost","Estimate"],2)
    results_sanc_pols$Competence[results_sanc_pols$Countries == "Pooled" & results_sanc_pols$Effect == "Conflict x Cost"] <- round(pool16$coefficients["territory:cost","Estimate"],2)
    
    results_sanc_pols$c_se <- NA
    results_sanc_pols$c_se[results_sanc_pols$Countries == "Spain-Portugal" & results_sanc_pols$Effect == "Conflict"] <- round(sp15$coefficients["territory","Std. Error"],2)
    results_sanc_pols$c_se[results_sanc_pols$Countries == "Spain-Portugal" & results_sanc_pols$Effect == "Cost"]     <- round(sp15$coefficients["cost","Std. Error"],2)
    results_sanc_pols$c_se[results_sanc_pols$Countries == "Brazil-Bolivia" & results_sanc_pols$Effect == "Conflict"] <- round(bb15$coefficients["territory","Std. Error"],2)
    results_sanc_pols$c_se[results_sanc_pols$Countries == "Brazil-Bolivia" & results_sanc_pols$Effect == "Cost"]     <- round(bb15$coefficients["cost","Std. Error"],2)
    
    results_sanc_pols$c_se[results_sanc_pols$Countries == "Pooled" & results_sanc_pols$Effect == "Conflict"] <- round(pool15$coefficients["territory","Std. Error"],2)
    results_sanc_pols$c_se[results_sanc_pols$Countries == "Pooled" & results_sanc_pols$Effect == "Cost"]     <- round(pool15$coefficients["cost","Std. Error"],2)
    results_sanc_pols$c_se[results_sanc_pols$Countries == "Spain-Portugal" & results_sanc_pols$Effect == "Conflict x Cost"] <- round(sp16$coefficients["territory:cost","Std. Error"],2)
    results_sanc_pols$c_se[results_sanc_pols$Countries == "Brazil-Bolivia" & results_sanc_pols$Effect == "Conflict x Cost"] <- round(bb16$coefficients["territory:cost","Std. Error"],2)
    results_sanc_pols$c_se[results_sanc_pols$Countries == "Pooled" & results_sanc_pols$Effect == "Conflict x Cost"] <- round(pool16$coefficients["territory:cost","Std. Error"],2)
    
    results_sanc_pols$lb_c <- results_sanc_pols$Competence - results_sanc_pols$c_se*1.616
    results_sanc_pols$ub_c <- results_sanc_pols$Competence + results_sanc_pols$c_se*1.616
    
    results_sanc_pols
    
    results_sanc_pols$Effect_c <- NA
    results_sanc_pols$Effect_c[results_sanc_pols$Effect == "Conflict"] <- "C"
    results_sanc_pols$Effect_c[results_sanc_pols$Effect == "Cost"] <- "B"
    results_sanc_pols$Effect_c[results_sanc_pols$Effect == "Conflict x Cost"] <- "A"
    
  # Combine Political Effects
    results_sanc_pols
    names(results_sanc_pols)
    
    results_h <- results_sanc_pols |>
      dplyr::select(Countries, Effect, Handling, h_se, lb_h, ub_h, Effect_h) |>
      rename(beta = Handling, se = h_se, lb = lb_h, ub = ub_h, order = Effect_h)
    
    results_h$Outcome <- "Handling"
    
    results_a <- results_sanc_pols |>
      dplyr::select(Countries, Effect, Approval, a_se, lb_a, ub_a, Effect_a) |>
      rename(beta = Approval, se = a_se, lb = lb_a, ub = ub_a, order = Effect_a)
    
    results_a$Outcome <- "Approval"
    
    results_c <- results_sanc_pols |>
      dplyr::select(Countries, Effect, Competence, c_se, lb_c, ub_c, Effect_c) |>
      rename(beta = Competence, se = c_se, lb = lb_c, ub = ub_c, order = Effect_c)
    
    results_c$Outcome <- "Competence"
    
    results_c
    
    pol_fx <- rbind(results_h,
                    results_a,
                    results_c)
    
    pdf(file="politics.pdf", height=10, width=16)
    
    ggplot(pol_fx, aes(x = beta, y = order, color = Countries)) + 
      theme_bw() +
      facet_wrap(~ Outcome) +
      geom_point(size = 5, position = position_dodge2(w = 0.75), shape = 15) +
      geom_linerange(aes(xmin = lb, xmax = ub), 
                     position = position_dodge2(w = 0.75),
                     lwd = 1.5) +
      scale_color_manual(values=c("lightgray", "darkgray", "black")) +
      labs(x = "\n Marginal Effect of Experimental Factors") +
      scale_y_discrete(name = "Effects \n",
                       labels = c("C" = "   Precipitant = \n Invasion", "B" = " Costly \n Sanctions", "A" = "Precipitant x Cost")) +
      theme(axis.title = element_text(size = 20)) +
      theme(axis.text = element_text(size = 20)) +
      theme(legend.title = element_text(size = 20, face = "bold"),
            legend.text = element_text(size = 20),
            strip.text = element_text(size = 20),
            legend.title.align=0.5) +
      theme(axis.text.y = element_text(angle=90, hjust = .5)) +
      geom_vline(xintercept = 0, lwd = 1.2, lty = 2)
    
    dev.off()
    
# Appendix: Results -------------------------------------------------------

  # Bolivia-Brazil (Support)
    bb1 <- summary(glm(sanctionynn ~ territory + cost, 
                       family = "binomial", data = bb))
    bb2 <- summary(glm(sanctionynn ~ territory + cost + territory:cost, 
                       family = "binomial", data = bb))
    bb3 <- summary(glm(sanctionynn ~ territory + cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       family = "binomial", data = bb))
    bb4 <- summary(glm(sanctionynn ~ territory + cost + territory:cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       family = "binomial", data = bb))
    
    bb1mod <- glm(sanctionynn ~ territory + cost, 
                  family = "binomial", data = bb)
    bb2mod <- glm(sanctionynn ~ territory + cost + territory:cost, 
                  family = "binomial", data = bb)
    bb3mod <- glm(sanctionynn ~ territory + cost +
                    as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                  family = "binomial", data = bb)
    bb4mod <- glm(sanctionynn ~ territory + cost + territory:cost +
                    as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                  family = "binomial", data = bb)
    
    texreg::texreg(list(bb1mod,bb2mod,bb3mod,bb4mod),
                   stars = c(.01,.05,.1),
                   digits = 3,
                   custom.coef.map = list("territory" = "Precipitant",
                                          "cost" = "Cost",
                                          "as.factor(partisan3)1" = "Independent",
                                          "as.factor(partisan3)2" = "Republican",
                                          "k40" = "Income",
                                          "millenial" = "Millenial",
                                          "white" = "White",
                                          "college" = "College Educated",
                                          "hinf" = "High Information",
                                          "salary" = "Employed for Salary",
                                          "territory:cost" = "Precipitant x Cost",
                                          "(Intercept)" = "Constant"))  
    
  # Bolivia-Brazil (Handling)
    bb5 <- summary(lm(handlingn ~ territory + cost, data = bb))
    bb6 <- summary(lm(handlingn ~ territory + cost + territory:cost, data = bb))
    bb7 <- summary(lm(handlingn ~ territory + cost +
                        as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                      data = bb))
    bb8 <- summary(lm(handlingn ~ territory + cost + territory:cost +
                        as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                      data = bb))
    
    bb5mod <- lm(handlingn ~ territory + cost, data = bb)
    bb6mod <- lm(handlingn ~ territory + cost + territory:cost, data = bb)
    bb7mod <- lm(handlingn ~ territory + cost +
                   as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                 data = bb)
    bb8mod <- lm(handlingn ~ territory + cost + territory:cost +
                   as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                 data = bb)
    
    texreg::texreg(list(bb5mod,bb6mod,bb7mod,bb8mod),
                   stars = c(.01,.05,.1),
                   digits = 3,
                   custom.coef.map = list("territory" = "Precipitant",
                                          "cost" = "Cost",
                                          "as.factor(partisan3)1" = "Independent",
                                          "as.factor(partisan3)2" = "Republican",
                                          "k40" = "Income",
                                          "millenial" = "Millenial",
                                          "white" = "White",
                                          "college" = "College Educated",
                                          "hinf" = "High Information",
                                          "salary" = "Employed for Salary",
                                          "territory:cost" = "Precipitant x Cost",
                                          "(Intercept)" = "Constant"))  
    
  # Bolivia-Brazil (Approval)  
    bb9  <- summary(lm(approvaln ~ territory + cost, data = bb))
    bb10 <- summary(lm(approvaln ~ territory + cost + territory:cost, data = bb))
    bb11 <- summary(lm(approvaln ~ territory + cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       data = bb))
    bb12 <- summary(lm(approvaln ~ territory + cost + territory:cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       data = bb))
    
    bb9mod  <- lm(approvaln ~ territory + cost, data = bb)
    bb10mod <- lm(approvaln ~ territory + cost + territory:cost, data = bb)
    bb11mod <- lm(approvaln ~ territory + cost +
                    as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                  data = bb)
    bb12mod <- lm(approvaln ~ territory + cost + territory:cost +
                    as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                  data = bb)
    
    texreg::texreg(list(bb9mod,bb10mod,bb11mod,bb12mod),
                   stars = c(.01,.05,.1),
                   digits = 3,
                   custom.coef.map = list("territory" = "Precipitant",
                                          "cost" = "Cost",
                                          "as.factor(partisan3)1" = "Independent",
                                          "as.factor(partisan3)2" = "Republican",
                                          "k40" = "Income",
                                          "millenial" = "Millenial",
                                          "white" = "White",
                                          "college" = "College Educated",
                                          "hinf" = "High Information",
                                          "salary" = "Employed for Salary",
                                          "territory:cost" = "Precipitant x Cost",
                                          "(Intercept)" = "Constant"))  
    
  # Bolivia-Brazil (Competence)
    bb13 <- summary(lm(competencen ~ territory + cost, data = bb))
    bb14 <- summary(lm(competencen ~ territory + cost + territory:cost, data = bb))
    bb15 <- summary(lm(competencen ~ territory + cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       data = bb))
    bb16 <- summary(lm(competencen ~ territory + cost + territory:cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       data = bb))
    
    bb13mod <- summary(lm(competencen ~ territory + cost, data = bb))
    bb14mod <- summary(lm(competencen ~ territory + cost + territory:cost, data = bb))
    bb15mod <- summary(lm(competencen ~ territory + cost +
                            as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                          data = bb))
    bb16mod <- summary(lm(competencen ~ territory + cost + territory:cost +
                            as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                          data = bb))
    
    texreg::texreg(list(bb13mod,bb14mod,bb15mod,bb16mod),
                   stars = c(.01,.05,.1),
                   digits = 3,
                   custom.coef.map = list("territory" = "Precipitant",
                                          "cost" = "Cost",
                                          "as.factor(partisan3)1" = "Independent",
                                          "as.factor(partisan3)2" = "Republican",
                                          "k40" = "Income",
                                          "millenial" = "Millenial",
                                          "white" = "White",
                                          "college" = "College Educated",
                                          "hinf" = "High Information",
                                          "salary" = "Employed for Salary",
                                          "territory:cost" = "Precipitant x Cost",
                                          "(Intercept)" = "Constant")) 
    
    
  # Spain-Portugal (Support)
    sp1 <- summary(glm(sanctionynn ~ territory + cost, 
                       family = "binomial", data = sp))
    sp2 <- summary(glm(sanctionynn ~ territory + cost + territory:cost, 
                       family = "binomial", data = sp))
    sp3 <- summary(glm(sanctionynn ~ territory + cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       family = "binomial", data = sp))
    sp4 <- summary(glm(sanctionynn ~ territory + cost + territory:cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       family = "binomial", data = sp))
    
    sp1mod <- glm(sanctionynn ~ territory + cost, 
                  family = "binomial", data = sp)
    sp2mod <- glm(sanctionynn ~ territory + cost + territory:cost, 
                  family = "binomial", data = sp)
    sp3mod <- glm(sanctionynn ~ territory + cost +
                    as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                  family = "binomial", data = sp)
    sp4mod <- glm(sanctionynn ~ territory + cost + territory:cost +
                    as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                  family = "binomial", data = sp)
    
    
    texreg::texreg(list(sp1mod,sp2mod,sp3mod,sp4mod),
                   stars = c(.01,.05,.1),
                   digits = 3,
                   custom.coef.map = list("territory" = "Precipitant",
                                          "cost" = "Cost",
                                          "as.factor(partisan3)1" = "Independent",
                                          "as.factor(partisan3)2" = "Republican",
                                          "k40" = "Income",
                                          "millenial" = "Millenial",
                                          "white" = "White",
                                          "college" = "College Educated",
                                          "hinf" = "High Information",
                                          "salary" = "Employed for Salary",
                                          "territory:cost" = "Precipitant x Cost",
                                          "(Intercept)" = "Constant"))  
    
  # Spain-Portugal (Handling)
    sp5 <- summary(lm(handlingn ~ territory + cost, data = sp))
    sp6 <- summary(lm(handlingn ~ territory + cost + territory:cost, data = sp))
    sp7 <- summary(lm(handlingn ~ territory + cost +
                        as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                      data = sp))
    sp8 <- summary(lm(handlingn ~ territory + cost + territory:cost +
                        as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                      data = sp))
    
    sp5mod <- lm(handlingn ~ territory + cost, data = sp)
    sp6mod <- lm(handlingn ~ territory + cost + territory:cost, data = sp)
    sp7mod <- lm(handlingn ~ territory + cost +
                   as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                 data = sp)
    sp8mod <- lm(handlingn ~ territory + cost + territory:cost +
                   as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                 data = sp)
    
    texreg::texreg(list(sp5mod,sp6mod,sp7mod,sp8mod),
                   stars = c(.01,.05,.1),
                   digits = 3,
                   custom.coef.map = list("territory" = "Precipitant",
                                          "cost" = "Cost",
                                          "as.factor(partisan3)1" = "Independent",
                                          "as.factor(partisan3)2" = "Republican",
                                          "k40" = "Income",
                                          "millenial" = "Millenial",
                                          "white" = "White",
                                          "college" = "College Educated",
                                          "hinf" = "High Information",
                                          "salary" = "Employed for Salary",
                                          "territory:cost" = "Precipitant x Cost",
                                          "(Intercept)" = "Constant"))  
    
  # Spain-Portugal (Approval)  
    sp9  <- summary(lm(approvaln ~ territory + cost, data = sp))
    sp10 <- summary(lm(approvaln ~ territory + cost + territory:cost, data = sp))
    sp11 <- summary(lm(approvaln ~ territory + cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       data = sp))
    sp12 <- summary(lm(approvaln ~ territory + cost + territory:cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       data = sp))
    
    sp9mod  <-lm(approvaln ~ territory + cost, data = sp)
    sp10mod <-lm(approvaln ~ territory + cost + territory:cost, data = sp)
    sp11mod <-lm(approvaln ~ territory + cost +
                   as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                 data = sp)
    sp12mod <-lm(approvaln ~ territory + cost + territory:cost +
                   as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                 data = sp)
    
    texreg::texreg(list(sp9mod,sp10mod,sp11mod,sp12mod),
                   stars = c(.01,.05,.1),
                   digits = 3,
                   custom.coef.map = list("territory" = "Precipitant",
                                          "cost" = "Cost",
                                          "as.factor(partisan3)1" = "Independent",
                                          "as.factor(partisan3)2" = "Republican",
                                          "k40" = "Income",
                                          "millenial" = "Millenial",
                                          "white" = "White",
                                          "college" = "College Educated",
                                          "hinf" = "High Information",
                                          "salary" = "Employed for Salary",
                                          "territory:cost" = "Precipitant x Cost",
                                          "(Intercept)" = "Constant"))  
    
  # Spain-Portugal (Competence)  
    sp13 <- summary(lm(competencen ~ territory + cost, data = sp))
    sp14 <- summary(lm(competencen ~ territory + cost + territory:cost, data = sp))
    sp15 <- summary(lm(competencen ~ territory + cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       data = sp))
    sp16 <- summary(lm(competencen ~ territory + cost + territory:cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       data = sp))
    
    sp13mod <- lm(competencen ~ territory + cost, data = sp)
    sp14mod <- lm(competencen ~ territory + cost + territory:cost, data = sp)
    sp15mod <- lm(competencen ~ territory + cost +
                    as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                  data = sp)
    sp16mod <- lm(competencen ~ territory + cost + territory:cost +
                    as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                  data = sp)
    
    texreg::texreg(list(sp13mod,sp14mod,sp15mod,sp16mod),
                   stars = c(.01,.05,.1),
                   digits = 3,
                   custom.coef.map = list("territory" = "Precipitant",
                                          "cost" = "Cost",
                                          "as.factor(partisan3)1" = "Independent",
                                          "as.factor(partisan3)2" = "Republican",
                                          "k40" = "Income",
                                          "millenial" = "Millenial",
                                          "white" = "White",
                                          "college" = "College Educated",
                                          "hinf" = "High Information",
                                          "salary" = "Employed for Salary",
                                          "territory:cost" = "Precipitant x Cost",
                                          "(Intercept)" = "Constant"))  
    
  # Pooled (Support)
    pool1 <- summary(glm(sanctionynn ~ territory + cost, 
                         family = "binomial", data = spbbdat))
    pool2 <- summary(glm(sanctionynn ~ territory + cost + territory:cost, 
                         family = "binomial", data = spbbdat))
    pool3 <- summary(glm(sanctionynn ~ territory + cost +
                           as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                         family = "binomial", data = spbbdat))
    pool4 <- summary(glm(sanctionynn ~ territory + cost + territory:cost +
                           as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                         family = "binomial", data = spbbdat))
    
    pool1mod <- glm(sanctionynn ~ territory + cost, 
                    family = "binomial", data = spbbdat)
    pool2mod <- glm(sanctionynn ~ territory + cost + territory:cost, 
                    family = "binomial", data = spbbdat)
    pool3mod <- glm(sanctionynn ~ territory + cost +
                      as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                    family = "binomial", data = spbbdat)
    pool4mod <- glm(sanctionynn ~ territory + cost + territory:cost +
                      as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                    family = "binomial", data = spbbdat)
    
    
    texreg::texreg(list(pool1mod,pool2mod,pool3mod,pool4mod),
                   stars = c(.01,.05,.1),
                   digits = 3,
                   custom.coef.map = list("territory" = "Precipitant",
                                          "cost" = "Cost",
                                          "as.factor(partisan3)1" = "Independent",
                                          "as.factor(partisan3)2" = "Republican",
                                          "k40" = "Income",
                                          "millenial" = "Millenial",
                                          "white" = "White",
                                          "college" = "College Educated",
                                          "hinf" = "High Information",
                                          "salary" = "Employed for Salary",
                                          "territory:cost" = "Precipitant x Cost",
                                          "(Intercept)" = "Constant")) 
    
  # Pooled (Handling)  
    pool5 <- summary(lm(handlingn ~ territory + cost, data = spbbdat))
    pool6 <- summary(lm(handlingn ~ territory + cost + territory:cost, data = spbbdat))
    pool7 <- summary(lm(handlingn ~ territory + cost +
                          as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                        data = spbbdat))
    pool8 <- summary(lm(handlingn ~ territory + cost + territory:cost +
                          as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                        data = spbbdat))
    
    
    pool5mod <- lm(handlingn ~ territory + cost, data = spbbdat)
    pool6mod <- lm(handlingn ~ territory + cost + territory:cost, data = spbbdat)
    pool7mod <- lm(handlingn ~ territory + cost +
                     as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                   data = spbbdat)
    pool8mod <- lm(handlingn ~ territory + cost + territory:cost +
                     as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                   data = spbbdat)
    
    texreg::texreg(list(pool5mod,pool6mod,pool7mod,pool8mod),
                   stars = c(.01,.05,.1),
                   digits = 3,
                   custom.coef.map = list("territory" = "Precipitant",
                                          "cost" = "Cost",
                                          "as.factor(partisan3)1" = "Independent",
                                          "as.factor(partisan3)2" = "Republican",
                                          "k40" = "Income",
                                          "millenial" = "Millenial",
                                          "white" = "White",
                                          "college" = "College Educated",
                                          "hinf" = "High Information",
                                          "salary" = "Employed for Salary",
                                          "territory:cost" = "Precipitant x Cost",
                                          "(Intercept)" = "Constant")) 
    
  # Pooled (Presidential Approval)  
    pool9  <- summary(lm(approvaln ~ territory + cost, data = spbbdat))
    pool10 <- summary(lm(approvaln ~ territory + cost + territory:cost, data = spbbdat))
    pool11 <- summary(lm(approvaln ~ territory + cost +
                           as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                         data = spbbdat))
    pool12 <- summary(lm(approvaln ~ territory + cost + territory:cost +
                           as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                         data = spbbdat))
    
    
    pool9mod  <- summary(lm(approvaln ~ territory + cost, data = spbbdat))
    pool10mod <- summary(lm(approvaln ~ territory + cost + territory:cost, data = spbbdat))
    pool11mod <- summary(lm(approvaln ~ territory + cost +
                              as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                            data = spbbdat))
    pool12mod <- summary(lm(approvaln ~ territory + cost + territory:cost +
                              as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                            data = spbbdat))
    
    
    texreg::texreg(list(pool9mod,pool10mod,pool11mod,pool12mod),
                   stars = c(.01,.05,.1),
                   digits = 3,
                   custom.coef.map = list("territory" = "Precipitant",
                                          "cost" = "Cost",
                                          "as.factor(partisan3)1" = "Independent",
                                          "as.factor(partisan3)2" = "Republican",
                                          "k40" = "Income",
                                          "millenial" = "Millenial",
                                          "white" = "White",
                                          "college" = "College Educated",
                                          "hinf" = "High Information",
                                          "salary" = "Employed for Salary",
                                          "territory:cost" = "Precipitant x Cost",
                                          "(Intercept)" = "Constant")) 
    
  # Pooled (Competence)  
    pool13 <- summary(lm(competencen ~ territory + cost, data = spbbdat))
    pool14 <- summary(lm(competencen ~ territory + cost + territory:cost, data = spbbdat))
    pool15 <- summary(lm(competencen ~ territory + cost +
                           as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                         data = spbbdat))
    pool16 <- summary(lm(competencen ~ territory + cost + territory:cost +
                           as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                         data = spbbdat))
    
    pool13mod <- lm(competencen ~ territory + cost, data = spbbdat)
    pool14mod <- lm(competencen ~ territory + cost + territory:cost, data = spbbdat)
    pool15mod <- lm(competencen ~ territory + cost +
                      as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                    data = spbbdat)
    pool16mod <- lm(competencen ~ territory + cost + territory:cost +
                      as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                    data = spbbdat)
    
    texreg::texreg(list(pool13mod,pool14mod,pool15mod,pool16mod),
                   stars = c(.01,.05,.1),
                   digits = 3,
                   custom.coef.map = list("territory" = "Precipitant",
                                          "cost" = "Cost",
                                          "as.factor(partisan3)1" = "Independent",
                                          "as.factor(partisan3)2" = "Republican",
                                          "k40" = "Income",
                                          "millenial" = "Millenial",
                                          "white" = "White",
                                          "college" = "College Educated",
                                          "hinf" = "High Information",
                                          "salary" = "Employed for Salary",
                                          "territory:cost" = "Precipitant x Cost",
                                          "(Intercept)" = "Constant")) 
    
# Auxiliary Analysis – Restricted Sample ----------------------------------

  # Restricted Data
    plot(density(spbbdat$duration,na.rm = TRUE))
    
    quantile(spbbdat$duration, probs = seq(0,1,.1), na.rm = TRUE)/60
    
    subset_spbbdat <- spbbdat  |> 
      filter(duration > 314 & duration < 636)
    
    subset_spdat <- sp  |> 
      filter(duration > 314 & duration < 636)
    
    subset_bbdat <- bb  |> 
      filter(duration > 314 & duration < 636)
    
  # % Change in Support
    sp3 <- summary(glm(sanctionynn ~ territory + cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       family = "binomial", data =subset_spdat))
    
    bb3 <- summary(glm(sanctionynn ~ territory + cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       family = "binomial", data =subset_bbdat))
    
    pool3 <- summary(glm(sanctionynn ~ territory + cost +
                           as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                         family = "binomial", data = subset_spbbdat))
    
    sp4 <- summary(glm(sanctionynn ~ territory + cost + territory:cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       family = "binomial", data =subset_spdat))
    
    pool4 <- summary(glm(sanctionynn ~ territory + cost + territory:cost +
                           as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                         family = "binomial", data = subset_spbbdat))
    
    bb4 <- summary(glm(sanctionynn ~ territory + cost + territory:cost +
                         as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                       family = "binomial", data =subset_bbdat))
    
    
    # Set Up Results Matrix  
    Countries <- c("Spain-Portugal","Brazil-Bolivia","Pooled")
    Effect <- c("Conflict","Cost","Conflict x Cost")
    
    results_sanc <- as_tibble(expand.grid("Countries" = Countries,
                                          "Effect" = Effect))
    results_sanc
    
    # Harvest Effects
    results_sanc$`Odds Ratios` <- NA
    
    results_sanc$`Odds Ratios`[results_sanc$Countries == "Spain-Portugal" & results_sanc$Effect == "Conflict"] <- round((exp(sp3$coefficients["territory","Estimate"]) - 1)*100,1)
    results_sanc$`Odds Ratios`[results_sanc$Countries == "Spain-Portugal" & results_sanc$Effect == "Cost"]     <- round((exp(sp3$coefficients["cost","Estimate"]) - 1)*100,1)
    
    results_sanc$`Odds Ratios`[results_sanc$Countries == "Brazil-Bolivia" & results_sanc$Effect == "Conflict"] <- round((exp(bb3$coefficients["territory","Estimate"]) - 1)*100,1)
    results_sanc$`Odds Ratios`[results_sanc$Countries == "Brazil-Bolivia" & results_sanc$Effect == "Cost"]     <- round((exp(bb3$coefficients["cost","Estimate"]) - 1)*100,1)
    
    results_sanc$`Odds Ratios`[results_sanc$Countries == "Pooled" & results_sanc$Effect == "Conflict"] <- round((exp(pool3$coefficients["territory","Estimate"]) - 1)*100,1)
    results_sanc$`Odds Ratios`[results_sanc$Countries == "Pooled" & results_sanc$Effect == "Cost"]     <- round((exp(pool3$coefficients["cost","Estimate"]) - 1)*100,1)
    
    results_sanc$`Odds Ratios`[results_sanc$Countries == "Spain-Portugal" & results_sanc$Effect == "Conflict x Cost"] <- round((exp(sp4$coefficients["territory:cost","Estimate"]) - 1)*100,1)
    results_sanc$`Odds Ratios`[results_sanc$Countries == "Brazil-Bolivia" & results_sanc$Effect == "Conflict x Cost"] <- round((exp(bb4$coefficients["territory:cost","Estimate"]) - 1)*100,1)
    results_sanc$`Odds Ratios`[results_sanc$Countries == "Pooled" & results_sanc$Effect == "Conflict x Cost"] <- round((exp(pool4$coefficients["territory:cost","Estimate"]) - 1)*100,1)
    
    results_sanc$se <- NA
    results_sanc$se[results_sanc$Countries == "Spain-Portugal" & results_sanc$Effect == "Conflict"] <- round((exp(sp3$coefficients["territory","Std. Error"]) - 1)*100,1)
    results_sanc$se[results_sanc$Countries == "Spain-Portugal" & results_sanc$Effect == "Cost"]     <- round((exp(sp3$coefficients["cost","Std. Error"]) - 1)*100,1)
    results_sanc$se[results_sanc$Countries == "Brazil-Bolivia" & results_sanc$Effect == "Conflict"] <- round((exp(bb3$coefficients["territory","Std. Error"]) - 1)*100,1)
    results_sanc$se[results_sanc$Countries == "Brazil-Bolivia" & results_sanc$Effect == "Cost"]     <- round((exp(bb3$coefficients["cost","Std. Error"]) - 1)*100,1)
    
    results_sanc$se[results_sanc$Countries == "Pooled" & results_sanc$Effect == "Conflict"] <- round((exp(pool3$coefficients["territory","Std. Error"]) - 1)*100,1)
    results_sanc$se[results_sanc$Countries == "Pooled" & results_sanc$Effect == "Cost"]     <- round((exp(pool3$coefficients["cost","Std. Error"]) - 1)*100,1)
    results_sanc$se[results_sanc$Countries == "Spain-Portugal" & results_sanc$Effect == "Conflict x Cost"] <- round((exp(sp4$coefficients["territory:cost","Std. Error"]) - 1)*100,1)
    results_sanc$se[results_sanc$Countries == "Brazil-Bolivia" & results_sanc$Effect == "Conflict x Cost"] <- round((exp(bb4$coefficients["territory:cost","Std. Error"]) - 1)*100,1)
    results_sanc$se[results_sanc$Countries == "Pooled" & results_sanc$Effect == "Conflict x Cost"] <- round((exp(pool4$coefficients["territory:cost","Std. Error"]) - 1)*100,1)
    
    results_sanc$lb <- NA
    
    results_sanc$lb <- results_sanc$`Odds Ratios` - results_sanc$se*1.68
    results_sanc$ub <- results_sanc$`Odds Ratios` + results_sanc$se*1.68
    
    results_sanc$Effect_n <- NA
    results_sanc$Effect_n[results_sanc$Effect == "Conflict"] <- "C"
    results_sanc$Effect_n[results_sanc$Effect == "Cost"] <- "B"
    results_sanc$Effect_n[results_sanc$Effect == "Conflict x Cost"] <- "A"
    
    results_sanc
    
    #setwd("../Plots/")
    
    pdf(file="subeffectsrr.pdf", height=8, width=20)
    
    ggplot(results_sanc, aes(x = `Odds Ratios`, y = Effect_n, color = Countries)) + 
      theme_bw() +
      geom_point(size = 5, position = position_dodge2(w = 0.75), shape = 15) +
      geom_linerange(aes(xmin = lb, xmax = ub), 
                     position = position_dodge2(w = 0.75),
                     lwd = 1.5) +
      scale_color_manual(values=c("lightgray", "darkgray", "black")) +
      labs(x = "\n % Change in Baseline Support for Sanctions") +
      scale_y_discrete(name = "Hypothesized Effects \n",
                       labels = c("C" = "Precipitant =\n Invasion", "B" = "Costly \n Sanctions", "A" = "Precipitant = Invasion \n x Costly Sanctions")) +
      theme(axis.title = element_text(size = 24)) +
      theme(axis.text = element_text(size = 20)) +
      theme(legend.title = element_text(size = 20, face = "bold"),
            legend.text = element_text(size = 20),
            legend.title.align=0.5) +
      theme(axis.text.y = element_text(angle=90, hjust = .5)) +
      geom_vline(xintercept = 0, lwd = 1.2, lty = 2) +
      xlim(-275,275)
    
    dev.off()
    
  # Predicted Probabilities  
    summary(pool4mod_sub <- glm(sanctionynn ~ territory + cost + 
                                  as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                                family = "binomial", data = subset_spbbdat))
    
    summary(sp4mod_sub <- glm(sanctionynn ~ territory + cost + 
                                as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                              family = "binomial", data = subset_spdat))
    
    summary(bb4mod_sub <- glm(sanctionynn ~ territory + cost + 
                                as.factor(partisan3) + k40 + millenial + white + college + hinf + salary, 
                              family = "binomial", data = subset_bbdat))
    
    mean(spbbdat$sanctionynn, na.rm = T)
    
    new_data_pool <- expand.grid(
      territory = c(0, 1),
      cost = c(0, 1),
      partisan3 = c(1),  # Assuming the first level as default
      k40 = mean(spbbdat$k40, na.rm = TRUE),     # Replace NULL with the mean of k40
      millenial = mean(spbbdat$millenial, na.rm = TRUE),
      white = mean(spbbdat$white, na.rm = TRUE),
      college = mean(spbbdat$college, na.rm = TRUE),
      hinf = mean(spbbdat$hinf, na.rm = TRUE),
      salary = mean(spbbdat$salary, na.rm = TRUE)
    )
    
    new_data_sp <- expand.grid(
      territory = c(0, 1),
      cost = c(0, 1),
      partisan3 = c(1),  # Assuming the first level as default
      k40 = mean(spbbdat$k40, na.rm = TRUE),     # Replace NULL with the mean of k40
      millenial = mean(spbbdat$millenial, na.rm = TRUE),
      white = mean(spbbdat$white, na.rm = TRUE),
      college = mean(spbbdat$college, na.rm = TRUE),
      hinf = mean(spbbdat$hinf, na.rm = TRUE),
      salary = mean(spbbdat$salary, na.rm = TRUE)
    )
    
    new_data_bb <- expand.grid(
      territory = c(0, 1),
      cost = c(0, 1),
      partisan3 = c(1),  # Assuming the first level as default
      k40 = mean(spbbdat$k40, na.rm = TRUE),     # Replace NULL with the mean of k40
      millenial = mean(spbbdat$millenial, na.rm = TRUE),
      white = mean(spbbdat$white, na.rm = TRUE),
      college = mean(spbbdat$college, na.rm = TRUE),
      hinf = mean(spbbdat$hinf, na.rm = TRUE),
      salary = mean(spbbdat$salary, na.rm = TRUE)
    )
    
    # Predict probabilities
    predicted_probabilities_pooled <- predict(pool4mod_sub, newdata = new_data_pool, type = "response", se.fit = TRUE)
    predicted_probabilities_sp <- predict(sp4mod_sub, newdata = new_data_pool, type = "response", se.fit = TRUE)
    predicted_probabilities_bb <- predict(bb4mod_sub, newdata = new_data_pool, type = "response", se.fit = TRUE)
    
    # Add probabilities to the new_data dataframe
    new_data_pool$pred <- predicted_probabilities_pooled$fit
    new_data_pool$se <- predicted_probabilities_pooled$se.fit
    new_data_pool$lb <- new_data_pool$pred - 1.645*new_data_pool$se
    new_data_pool$ub <- new_data_pool$pred + 1.645*new_data_pool$se
    
    new_data_pool$Countries <- NA
    new_data_pool$Countries <- "Pooled"
    
    new_data_pool
    
    new_data_sp$pred <- predicted_probabilities_sp$fit
    new_data_sp$se <- predicted_probabilities_sp$se.fit
    new_data_sp$lb <- new_data_sp$pred - 1.645*new_data_sp$se
    new_data_sp$ub <- new_data_sp$pred + 1.645*new_data_sp$se
    
    new_data_sp
    
    new_data_sp$Countries <- NA
    new_data_sp$Countries <- "Spain-Portugal"
    
    new_data_bb$pred <- predicted_probabilities_bb$fit
    new_data_bb$se <- predicted_probabilities_bb$se.fit
    new_data_bb$lb <- new_data_bb$pred - 1.645*new_data_bb$se
    new_data_bb$ub <- new_data_bb$pred + 1.645*new_data_bb$se
    
    new_data_bb$Countries <- NA
    new_data_bb$Countries <- "Brazil-Bolivia"
    
    new_data_bb
    
    pred_data <- rbind(new_data_pool,new_data_sp,new_data_bb)
    
    pred_data$Cost <- NA
    pred_data$Cost[pred_data$cost == 1] <- "Costly"
    pred_data$Cost[pred_data$cost == 0] <- "Not Costly"
    
    # Convert 'territory' to a factor and set levels
    pred_data$territory <- factor(pred_data$territory, levels = unique(pred_data$territory))
    
    pred_data$Precipitant <- NA
    
    pred_data$Precipitant[pred_data$territory == 0] <- "Trade"
    pred_data$Precipitant[pred_data$territory == 1] <- "Invasion"
    
    # Plot with the sorted data
    pdf(file="restrictedprprobsrr.pdf", height=8, width=20)
    
    ggplot(pred_data, aes(x = Precipitant, y = pred, color = Cost)) +
      theme_bw() +
      geom_errorbar(aes(ymin = lb, ymax = ub), width = 0.5, position = position_dodge(width = 0.6), lwd = 1) +
      geom_point(size = 7, position = position_dodge(width = 0.6)) +
      scale_color_manual(values=c("lightgray", "black")) +
      facet_wrap(~Countries) +
      theme(axis.title = element_text(size = 20)) +
      theme(axis.text = element_text(size = 20)) +
      theme(legend.title = element_text(size = 20, face = "bold"),
            legend.text = element_text(size = 20),
            strip.text = element_text(size = 20),
            legend.title.align=0.5) +
      theme(axis.text=element_text(size=18),
            axis.title=element_text(size=24),
            plot.margin = margin(0.5,0.5,1.1,1.1,"cm"),
            axis.title.x = element_text(vjust = -2),
            axis.title.y = element_text(vjust = 4)) +
      labs(x = "Precipitant", y = "Probability of Approving of Sanctions", color = " Sanctions \nCost") +
      geom_hline(yintercept = 0.5, lwd = 1.2, lty = 2) +
      ylim(0,1)
    
    dev.off()
    
    
    
